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FOREWORD 


This  report  gives  an  overview  of  the  orbit  computation  program 
CELEST  used  by  the  Naval  Surface  Weapons  Center  in  processing  Doppler 
data  to  position  both  satellites  and  surface  based  Doppler  receivers. 
The  Doppler  system  has  been  in  use  since  1962  by  the  U.S.  Navy  and  has 
been  used  at  the  Naval  Surface  Weapons  Center,  Dahlgren,  Virginia 
primarily  for  Defense  Department  related  geodetic  studies. 

This  CELEST  computer  program  uses  raw  Doppler  data  to  determine 
satellite  orbits.  It  provides  diagnostic  information  on  the  quality 
of  the  orbits.  The  basic  technique  employed  is  one  of  weighted  least 
squares  where  the  data  is  edited  and  weighted  within  the  program.  An 
iterative  capability  exists  for  nonlinear  problems.  Trajectories  are 
formed  by  directly  integrating  the  equations  of  motion  in  an  inertial 
frame.  The  force  equation  has  components  due  to  Earth,  Sun  and  Moon 
gravity,  solar  radiation,  thrust,  atmospheric  drag,  solar  and  lunar 
tidal  distortion.  A satellite  frequency  offset  error  can  be  determined 
and  the  program  has  the  facility  for  determining  unknown  receiver 
locations.  The  computer  program  occupies  130K  octal  units  of  memory, 
is  structured  as  nine  major  overlays  and  is  completely  written  in 
fortran.  The  program  is  primarily  operated  on  a sixty  bit  CDC  6700 
computer  at  Dahlgren,  two  IBM  1108  computers  at  the  Defense  Mapping 
Agency  Aerospace  Center  and  Topographic  Center,  on  a CDC  6700  at 
Cambridge  Research  laboratory  and  an  SEL  86  located  at  Naval  Space 
Surveillance  Center,  Dahlgren,  Virginia. 

This  report  was  prepared  by  James  W.  O'Toole  and  reviewed  by 
Richard  J.  Anderle  of  the  Astronautics  and  Geodesy  Division. 
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MEASUREMENTS 


The  majority  of  Doppler  measurements  today  are  taken  by  Geoceiver 
or  Geoceiver  type  equipment.  This  equipment  integrates  the  Doppler 
effect  on  the  transmitted  frequency  fs  over  approximately  30  seconds 
to  obtain  a measurement  equivalent  to  range  difference.  Letting  fr 
be  the  receiver  generated  frequency,  slightly  offset  from  fs,  and  p 
indicate  range  from  the  receiver  to  the  satellite  we  have 

(1)  Doppler  = Nc  = fr  - fs(l  - p/c) 

Integrating  fives 

(2)  p2  - pi  = c/fs  [Nc  - (fr  - fs)  (t2  - t-, ) ] 

where  t2  - t-|  = 30  sec. 

c = velocity  of  light 

p = | Remission  time)  - r(reception  ti me ) j 
satellite  receiver 

Older  Doppler  equipment  which  integrates  over  a time  period  less  than 
a second  treats  the  data  as  instantaneous  range  rate.  This  is  done 
by  solving  equation  (1)  for  p and  assigning  the  midpoint  of  the 
integration  interval  as  the  time  of  observation. 

STRUCTURAL  OVERVIEW 

The  basic  program  modules  are  indicated  in  Diagram  1.  Coordinates 
of  the  sun  and  moon  are  retained  at  one  day  and  one  half  day  intervals 
respectively  on  the  Sun-Moon  file.  They  are  in  the  Mean  Inertial  System 
of  1950.0  and  ephemeris  time.  In  addition  values  for  the  inertial  to 
earth  fixed  coordinate  transformation  are  retained  on  this  file.  These 
values  are  the  nutation  in  longitude  and  obliquity  of  the  ecliptic  of 
the  sun,  Bessel ian  day  numbers  and  the  equation  of  the  equinoxes.  The 
Satellite  Table  contains  offset  frequency  values  to  the  broadcast 
frequency  of  individual  satellites.  The  Gravity  file  contains 
spherical  harmonic  coefficients. 


r 


Sun-Moon 

Satellite  Table  p These  files  must  always  be  attached  to  the  program. 
Gravity  j 


Diagram  1 

Basic  Program  Modules 
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The  program  utilizes  various  pre-processors  as  each  one  is  designed 
to  process  a specific  type  of  data.  Pre-processors  perforin  the 
function  of  preliminary  data  editing,  time  correcting  and  placing  the 
data  in  a specific  format  on  the  Data  Prep  fi le. 

Final  data  point  editing  and  the  assignment  of  weights  take  place 
in  the  filter  module.  This  function  is  performed  by  fitting  a 
satellite  trajectory  to  one  pass  of  data  at  a time.  Points  are  re- 
moved which  have  a residual,  from  the  fitted  orbit,  greater  than  twice 
the  computed  standard  deviation.  This  procedure  is  iterated  until  no 
further  outlying  points  are  identified.  At  this  point  the  RMS  of 
residuals,  from  the  fitted  orbit,  is  assigned  as  the  standard  deviation 
of  an  observation.  The  reciprocal  of  this  value  is  assigned  as  a weight. 
Least  square  normal  equations  are  made  up  using  the  untagged  points  and 
their  weights. 

The  normal  equations  contain  a subset  of  six  osculating  orbital 
elements,  one  drag  scaling  constant,  three  thrust  parameters,  a 
radiation  parameter,  three  receiver  station  coordinate  parameters,  a 
refraction  correction  parameter,  frequency  bias  and  frequency  drift 
parameters.  The  subset  is  determined  by  the  nature  of  the  problem 
being  solved.  The  normal  equation,  the  time  of  closest  approach  and 
the  RWS  of  residuals  are  written  to  a file  (Pass  Matrix  file). 

Ordinarily  the  Pass  Matrix  file  would  be  the  only  thing  required 
to  complete  processing  as  the  remaining  task  is  that  of  summing  matrices 
and  performing  an  inversion.  The  solution  area  is  designed  however, 
to  determine  solutions  over  various  time  domains  with  an  option  for  drag 
segmentation.  This  process  requires  access  to  partial  derivatives 
from  the  trajectory.  Segmentation  is  a process  whereby  the  effect  of 
perturbing  x number  of  drag  segments  independently  is  desired  as  part 
of  the  solution.  The  program  accomplishes  this  by  storing  the  effect 
of  perturbing  any  number  of  drag  segments  an  equal  amount  in  the  pass 
normal  matrix  via 


3 Data = 3D_ 

3 Drag  Parameter  3Cq 


and  using  the  trajectory  partials  in  the  solution  area  to  transform  the 
pass  matrix  value 


aD 

3C'd 


into  its  associated  values 


3D 


Each  3D  represents  the 
3Coi 
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perturbation  in  data  due  to  independently  perturbing  the  drag 
coefficient  in  the  i th  segment  of  the  drag  force.  This  permits  a 
flexible  treatment  of  drag  and  is  especially  designed  for  short  arc 
processing. 

The  final  step  in  the  orbit  determination  process  is  generating  a 
refined  ephemeris.  This  is  done  by  using  trajectory  partials  to 
compute  the  effect  on  the  reference  orbit  due  to  parameter  changes 
indicated  in  the  solution.  Trajectory  values  beyond  the  reference 
trajectory  time  span  must  of  course  be  generated  via  integration. 

The  central  theme  on  which  both  the  Celest  program  and  an 
independent  Celest  Station  Position  program  operates  is  the  data 
bank  concept  of  Diagram  2. 

It  is  pass  normal  matrices  which  are  saved  in  the  data  bank.  The 
matrices  become  the  data  for  most  work  although  raw  data  is  saved  for 
more  fundamental  studies.  By  a continued  use  of  matrices,  trajectories 
and  the  principle  of  first  order  matrix  adjustment,  due  to  a change  in 
the  reference  parameter  values,  orbit  refinement  and  station  positioning 
can  be  carried  on  with  a minimum  of  data  reprocessing. 

DATA  FILTERING 

Point  filtering  proceeds  with  the  philosophy  that  model  error 
during  a single  pass  can  be  removed  by  adjusting  the  orbit  in  the 
along  track  and  radial  directions  together  with  removing  frequency  and 
tropospheric  refraction  error.  The  ionospheric  error  has  been  re- 
moved at  the  receiving  station  by  gathering  dual  frequency  data  and 
eliminating  the  first  order  ionospheric  effect. 

The  filter  process  is  implemented  by  forming  a pass  normal  matrix 
on  a reference  orbit. 


^orbit,  orbit 

^orbit,  bias 

[~ Aorbi  t 

^orbi t 

®bias,  bias 

Abias 

Ebias  _ 

The  orbit  section  always  contains  six  osculating  orbital  elements 
at  some  epoch  time.  The  bias  section  contains  receiver  station 
coordinates,  a frequency  bias  parameter  which  measures  the  deviation 
of  the  satellite  frequency  offset  from  an  assumed  value  and  a re- 
fraction correction  giving  the  percent  deviation  from  the  Hopfield 
tropospheric  model.  A transformation  $ is  formed  whi ch  takes  normal 
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Raw  Data 


Trajectory  for  some 
reference  time  span 


Existing  Pass  Matrix  file  Pass  Matrix  file  of  new  data 


Merge 


Updated  Pass  Matrix  file 


Diagram  2 

Pass  Matrix  Data  Bank  Concept 
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equations  presented  in  osculating  orbital  elements  and  reforms  them 
in  position  velocity  components  resolved  in  a special  local  frame  at 
TCA.  The  local  frame  is  defined  by  the  station-satellite  range  ( p ) 
and  along  track  ( p ) vectors  at  the  satellites  time  of  closest 
approach. 

R = [p,  p,  pxp] 

where  the  * indicates  unit  vectors. 

The  epoch  transformation  going  from  the  osculating  orbital  element 
time  to  TCA  is 

if  (TCA)  = 3 (x  > x)  (TCA) 

3e0  (epoch  elements) 

The  transformation  <t>  is  then 

f = Rif  R = [§  °] 

A solution  is  determined  and  the  RWS  of  residuals  computed  on  the 
adjusted  orbit.  Data  point  weights  are  recomputed  by  dividing  their 
present  value  by  the  adjusted  RWS  of  residuals.  Individual  data 
points  are  tagged  if  the  product  of  their  weight  and  associated 
adjusted  residual  is  greater  than  two  (2  sigma  filtering).  The 
process  iterates  reforming  the  normal  equation  at  each  iteration  with 
the  untagged  points  and  improved  weights.  All  points  are  examined  at 
each  iteration  and  the  process  terminates  when  the  same  points  are 
tagged  on  successive  iterations.  The  final  normal  equation  and  the  RWS 
of  unadjusted  residuals  using  untagged  points  is  stored  on  the  Pass 
Matrix  file.  Basic  technical  procedures  used  in  this  work  are  given 
below. 


STATION  COORDINATE  CONVERSION 


where 


rSo  = earth  fixed  station  coordinates 

hs  = Station  geodetic  height  above  a reference  ellipsoid 

= eccentricity  of  the  reference  ellipsoid 

1 


= [(2  - f)f] 


2 - 


EL 


EL 


f = flattening 
EL  = oblateness 
= 298.25 

ae  = semi-major  axis  of  the  reference  ellipsoid 
= 6378.145  km 


Qs  = station  geodetic  latitude 

X0  = longitude  from  Greenwich  meridian 

A = ae/(  i25in2Qs)^ 

RANGE  AND  RANGE  RATE  COMPUTATION 


(4)  Range  = p(tr)  = r(te)  - rs(tr) 

where  tp  = tr  - (p*p)!s/c  tr  = reception  time  and  te  is  determined  by 


iteration  starting  with  te  = tr. 
(5)  p = r(te)  - rs(tr)  - £* 

where  rs  = u xr<. 


Jr(te)  - rSi(tr)].  rjtfi) 
c(p*p)1'5 

r-  = inertial  station  location 
= (ABCD)*  rSo 
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CD  = Precession  Nutation  transformation 


Earths  mean  sideral  rotation  rate  determined  by  the  time 
lapse  between  successive  transits  of  the  mean  equinox. 


= .7292115855  E(-04)  >~adians/sec . 
TIME  OF  CLOSEST  APPROACH  (TCA)  COMPUTATION 


Define 


t i+1  = ti  - (p*p)ti/VP*P  + P*p)ti 


Set  t]  = value  from  the  observation  file  and  iterate  until 

I ti+l 


- t^|  < .05  sec.  or  maximum  iteration  count  is  reached. 


COMPUTATION  OF  ZENITH  ANGLE 

(6)  Zv  = tan’1  ([1-(^*US)2]  /p*Us) 


Us  = (ABCD)*  /s1 


REFRACTION  MODEL  (H0PF1ELD) 


s^  computed  from  (3) 


with  hs=0 


nT] 

= dry  term 

E = 

Nf2 

= wet  term 

H = 

Ro 

= 1 rs0  1 

P = 

R1 

= R0  + 40.1  + . 1 49T 

V 

r2 

= R0  + 12.0 

= 

Ni 

= [(.776)  E-04]P/T|< 

H = 

n2 

= [(.373)  E-02]E/T|<2 

P = 

T = 

T = 15  - 6. 5h_ 


3 


* 


rsa  te 1 1 i te 
(n-l) 


Station  (r*r-k?7 


, f f rsatel 1 i te 
02)  AR  - k k (n_i) 


where 


kk  = o cos  Zy  [(I-pp*)p  -u  xp]*Us 


Letting  RQ  = rstation  and  Ri  = rsateiiite  we  can  compute  range  and 
range  rate  corrections  by 


/Ki  c.  4 • 1 4 • 

(n-l)  r»dr  = ■r’ _Jkj V'_l_  (4)  (r*r-k2)J*’'(  k2-R?)  J 

R0(r*r-k2)'5  - R?>4  j,0  y+i 

Ri  2 a 

(Hi  ar  ,[ ("-.lU^r-  yJh y-L- (i-r-k^iJ-V-M)4-1 

r-k*)V2  ^(RR-R?)4  £,M-1  J 


DATA  TYPE  FORMULATION 


Letting  Dg  denote  range  difference  and  D7  denote  range  rate  we  have 
(15)  D 9 = [ [p  | - fb(t-TCA)  - — fb(t-TCA)2  + n+cR)AfR]t'1 


where  f00  is  an  input  quantity  usually  taken  to  be  E +(06)  so  that  the 
bias  solution  will  be  in  ppm. 

Partial  derivatives  are  given  by 


5 CD  -j  _ 3[]-j_i 

3q  3q  3q 


q = orbit  set  ( p) , 

station  set  (rrn), 
CR.  fb  or  *b- 


.H 


07) 

iLL  = 

p* 

ax_ 

3p 

3P 

(18) 

3[]  = 

-p* 

(ABCD)* 

8 rSo 

(19) 

ILL  = 

&fR 

3cr 

(20) 

oLL  = 

-c 

(t-TCA) 

3fb 

foo 

(21) 

iLL  = 

-c 

(t-TCA)2 

A 

The  vacuum  received  frequency  using  a transmitted  frequency  of  f,  is 
given  by  b 


(22)  fR  - f$ 

^ 1 +p  * r(te)  J 

where  P(tr)  = r(te)  - rs(tr) 

(te)  = tr-(P*PL  = emitted  time 
c 

tr  = received  time 

computing  (5  from  the  expression  for p and  taking  into  account  the 
change  in  te  as  a function  of  tr  gives 


h * rs(tr)  \ 
\ c ) 


(23)  p = ^~(te^ ^s( V)  + (p*/c)  y (V)r(te)  ~ (f5*/c)r(te)rs(tr) 


HP*r(te)  -/(p*r(te) 

c \ c 

1 + P*r(te) 

C 
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1 + p*  rs(tr)/c 
1 + P*  r (te)/c 

= II 


Thus  if  p is  computed  to  first  order  in  1/c  fr/fs  will  be  given  to 
second  order  in  1/c. 


1st 


(24)  P=rder  f’(te)  - fs(tr)  - (i5*/c)  (r(te)  - fs(tr))^(te) 

Adding  in  the  contribution  of  frequency  bias  and  refraction  gives  the 
formulation  for  range  rate  as 

(25)  D7  = fsO+fb/foo+fo0(t'TCA))  + 0+cR)AfR 


(25) 

°7  = 

This 

formula 

formula  (15) 

(26) 

3D7  = 
3p“ 

3D-, 

(27) 

7 = 

3rs 

bo 

(28) 

3D7 

3Cr 

(29) 

3D7  _ 
8fb 

(30) 

3D7 

A* . 


~h  [_^i(  i-«*) 


I p I 
p*< 


3r 

3p 


*—  ] 
3p 


Is  [ ( l-|5p*)+p*Xu)*]  (ABCD)* 

c I p I 

AfR 


307  = fs_  (l-p*p/c) 


1 00 


3fh 


(s  (l-p*p/c)  (t-TCA) 


'00 


The  primary  output  statistics  from  the  filtering  process  are  filtered 
noise  and  the  RWS  of  residuals. 


(31) 


RWS  = 


N, 


obs 


X wi  (Dobs-D)i/Nobs 
i=l 


1/2 
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r 


'I 


where  D is  computed  on  the  reference  ephemeris  and  u.  is  the  computed 
weight. 

Nobs 

r 2]  1/2 

y (Nobs^i ) | 


(32)  Filtered  Noise  = ^ 

i=l 


= standard  deviation  for  an  observation 
from  the  pass. 


As  the  RWS  of  adjusted  residuals  converges  to  one  during  the 
filtering  process  and  the  weights  are  approximately  constant  we  have 

^obs 

[ (Dobs  ” Da(jjusted^  /^obs^2  - ^ 

1 _ 


i = l 


or  RMS  of  adjusted  residuals  = — = standard  deviation 

w = Filtered  Noise 


SOLUTION  OF  NORMAL  EQUATIONS 


The  solution  area  solves  an  equation  containing  from  six  to  thirty 
nine  dynamic  parameters,  several  hundred  bias  parameters  and  several 
sets  of  station  coordinates.  To  accomplish  this  the  program  performs 
bias  and  station  parameter  elimination  in  order  to  keep  the  matrix 
requiring  inversion  under  dimension  of  forty.  Letting  "0"  stand  for 
the  orbit  (dynamic)  parameters  and  "b"  for  bias  we  have  the  normal 
equation 


— — 1 
Boo  B0b 

AFb 

B0~ 

_Q  * 
_Q 
CO 

O 

_Q 

CO 

1 

^Pb 

Jb_ 

The  bias  terms  are  eliminated  from  the  above  equation  and  the  elimination 
equations  are  saved  in  order  to  recover  bias  solutions. 


(33) 

(34) 


APb  = Bbb  (Eb-BboAPc>) 


(^oo'^ob^bb^bo^  APo  ^0  ^ob^bb^b 


REliminated 
or  B00 


ARo 


-El i mi nated 
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Just  prior  to  bias  elimination  pass  matrices  are  expanded  from 
their  basic  set  of  parameters  to  the  desired  solution  set.  This 
is  a major  task  in  the  solution  area.  The  solution  set  may  consist 
of  orbital  elements  at  a different  epoch,  multiple  drag  parameters 
due  to  drag  segmentation,  multiple  thrust,  polar  motion  parameters, 
up  to  twenty  gravity  parameters,  bias  parameters  and  station  coordinates. 

Subsequent  to  bias  elimination  pass  matrices  are  summed  to  form 
an  arc  normal  matrix.  Station  coordinate  sections  of  the  matrix  are 
summed  over  each  station  and  eliminated  just  prior  to  matrix  inversion. 
Station  solutions  are  obtained  by  backsubsti tution.  Direct  observation 
of  all  parameters  is  permitted  in  the  form  of  apriori  sigma  input  for 
each  parameter.  The  inverse  square  of  the  input  sigma  is  added  to  the 
diagonal  term  of  the  normal  matrix  prior  to  inversion.  This  corresponds 
to  an  observation  of  the  present  value  of  a parameter  with  the  input 
sigma  as  the  standard  error  of  observation. 

NAVIGATION  SOLUTIONS 

The  primary  diagnostic  output  is  a set  of  navigation  solutions. 

A navigation  solution  is  carried  out  for  each  pass  of  data  and  consists 
of  determining  the  receiver  motion,  in  the  along  track  and  radial 
directions,  which  minimizes  residuals  for  the  refined  ephemeris.  The 
receiver  motion  is  interpreted  as  satellite  ephemeris  error  as  the 
receiver  location  is  actually  well  known.  The  radial  and  along  track 
directions  used  are  the  range  and  range  rate  vectors  at  TCA  discussed 
earlier  under  filtering.  The  procedure  is  implemented  by  saving  parts 
of  the  bias  elimination  equations  (33) 

Bbb>  Bbo>  Eb 

and  adjusting  for  the  orbit  solution  &po  obtained  at  the  time  of  arc 
matrix  inversion. 


Bbb  Ab  = Eb  - Bbo  Apo 

The  Bbb  section  contains  earth  fixed  receiver  coordinates  which  are 
transformed  to  the  local  TCA  frame.  Solutions  are  generated  for  radial, 
along  track,  frequency  and  refraction  error.  Diagnostic  cards  con- 
taining this  information  and  filter  statistics  are  generated  at  this 
time. 
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STATION  DIAGNOSTIC  CARDS 


Station  Analysis  cards  can  be  punched  on  option  and  contain 
diagnostic  informaion  useful  to  the  satellite  tracking  stations. 

The  card  values  are  as  follows: 

1.  STA  - Station  number  extracted  from  the  header  message  of  the 
raw  data. 

2.  Time  (hr,  min)  - Time  of  the  first  data  point. 

3.  TCA  (sec)-  Time  of  closest  appraoch  of  the  satellite  obtained 
from  the  satellite  orbit  by  searching  for  the  point  in  time  where  the 
range  (p)  and  range  rate  (p)  vectors  are  perpendicular  to  each  other. 


4.  FREQ(mc)  - The  Q number  from  the  raw  data  header  message  is 
used  to  determine  a base  frequency.  The  frequency  given  is  this  base 
frequency  rounded  to  the  nearest  MHz. 

5.  EL  (deg)  - The  satellite  elevation  at  TCA  computed  from  the 
satel 1 i te  orbit. 

6.  PTS.  Good  - Total  number  of  points  left  after  passing 
through  the  Celest  point  filtering  process. 

a.  Points  are  filtered  out  in  the  Celest  Pre-Processor  if 
information  is  missing,  values  are  to  large  or  the  data  fails  a 
monotone  test. 

b.  Points  are  filtered  out  in  the  Celest  Filter  by  removing 
orbit  error  from  the  residuals  and  testing  against  2.0  sigma. 
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7.  Flit.  Noise  (cps/MHz  or  in)  - Filtered  noise  is  the  standard 
deviation  on  the  data  after  modeling  error  has  been  removed.  In  the 
case  of  range  rate  this  value  is  given  in  units  of  cps/MHz.  In  the 
cases  of  range  difference  the  value  is  given  in  units  of  meters. 


Before  modeling 
error  removed 


After  modeling 
error  removed 


Filtered  Noise  = RMS  of  residuals  after  modeling  error  is  removed. 

This  value  is  scaled  to  1 MHz  for  Doppler  data. 

8.  ELT.(m)  - This  is  the  along  track  navigation  error  determined 
from  the  refined  orbit.  Holding  the  orbit  fixed  the  station  is 
allowed  to  move  in  the  along  track  (p)  and  slant  range  (p)  directions, 
in  order  to  best  fit  the  data  from  the  pass.  Since  this  movement  is 
from  a known  position  the  result  is  tabulated  as  along  track  (ELT) 
and  slant  range  (ELR)  errors.  The  values  represent  a measurement  of 
how  well  the  final  refined  orbit  fits  the  data  of  a given  pass,  (see 
diagram  for  #3) 

9.  ELR.(m)  - Slant  range  error. 
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10.  DLT.  F.  (ppm)^)  - Delta  frequency  is  the  value  of  the 
frequency  bias  determined  in  the  above  navigation  solution.  Assuming 
that  the  satellite  frequency  has  a constant  bias  during  a given  pass, 
then  this  number  represents  that  bias  in  parts  per  million. 

11.  ACT  - Action  taken  in  the  course  of  point  filtering.  Action 
label  described  below. 

A - No  TCA 

B - Rejected  in  filter  because  to  many  points  were  filter  out. 

D - Rejected  on  TCA  zenith  angle  test. 

E - Pass  not  balanced,  i.e.  the  difference  between  the  number 
of  points  on  one  side  of  TCA  and  the  number  of  points  on 
the  other  is  greater  than  the  balanced  pass  tolerance. 

Passes  are  rejected  for  reasons  other  than  the  above,  in  the  Pre- 
Processor.  These  reasons  are  listed  in  Pre-Processor  under  Reject  Codes 
but  no  indication  is  given  on  the  Station  Analysis  cards. 

INTEGRATION 

For  ephemeris  computation  a variable  order  routine  is  used  on  a 
14  digit  machine.  Usually  the  order  is  set  to  ten  and  navigation 
type  satellites  use  a 60  sec  step  size.  The  ephemeris  can  be  computed 
to  a one  meter  accuracy  and  perturbations  of  orbit  constants  to  one  part 
in  106. 

The  integration  routine  is  a Gauss  Jackson  technique  using  backward 
differences  and  follows  the  basic  pattern  of  first  initializing  a 
backward  difference  table  to  the  order  (N)  of  the  process 


( ) A nominal  value  of  oscillator  frequency  offset  (in  parts  per 
million)  is  associated  with  each  satellite.  On  the  basis  of 
the  Doppler  data  from  a given  pass,  a correction,  DLT.  F.  is 
calculated.  The  corrected  absolute  offset,  in  parts  per 
million,  for  the  pass  is  the  algebraic  sum  of  the  nominal 
offset  and  DLT.  F. 

Absolute  offset  = Nominal  offset  + DLT.  F. 

= Avs  + DLT.  F. 

The  nominal  offset  is  always  negative.  For  example,  if  the 
nominal  offset  for  a satellite  is  -80  ppm,  a DLT.  F.  of  +.04 
would  indicate  an  absolute  offset  of  -80+.04  = -79.96  ppm. 


£ 
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then 


Extrapolates  the  difference  table  from  line  n to  n+1 . 

v-1  x(n)  = v"'*  x(n-l)  + 7°  x(n) 

V'2  x(n)  = x(n-l ) + v-lx(n) 

vN  x(n+l)  = vN  x(n) 

k v 

v x(ml)=  x(n ) + vK^x(n+l ) K=N-1 , . . .0 


2.  Compute j position  and  velocity  by 

N 


where 


Xn+l=h2v"2xn+h2  Y=0  cK 


2 


VK  y 

^ xn+ 1 


<n+l 


= h v 


-1-J 


N 


h 2 


K=0 


aK  ~ aKHl 

ao  " 1+al 

K+2 

n 

7^ 

II 

aeaK+2-e 

K 

aK- j 
J+j 

■k  - - 2 

J=1 

xn+l 

a0=l  a -|  =-l /2 


3.  Uses  the  force  function,  G,  from  x = G(x,x,t)  to  compute  x(n+l) 
and  determine  the  difference  between  the  computed  and  extrapolated 
values  of  acceleration.  The  backward  difference  table  is  adjusted  due 
to  this  difference. 


4.  The  process  described  in  2 and  3 is  continued  until  the  desired 
number  of  iterations  is  reached.  The  final  result  from  2 is  written 
on  the  trajectory  file  and  the  process  terminates  by  carrying  out  step  3. 

COORDINATE  AND  TIME  SYSTEMS 


The  reference  frame  used  for  integration  is  an  inertial  frame  defined 
by  the  mean  equator  and  equinox  of  1950.0  the  sun  and  moon  coordinates 
are  in  the  1950.0  system  using  Ephemeris  Time  (ET).  The  time  system 
for  Doppler  observations  is  Universal  Time  Coordinates  (UTC)  and  thus 
the  integration  time  is  UTC.  The  difference  between  UTC  and  ET  is 
presently  46.15  sec  and  is  not  presently  adjusted  for  when  using  the 
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and  moon  coordinates  in  the  force  model.  The  gravity  force  is  computed 
by  rotating  the  satellite  inertial  position  to  the  earth  fixed  frame 
aligned  with  the  Greenwich  meridian  and  using  the  instantaneous  earth 
spin  axis.  The  difference  between  UTC  and  UT1  is  taken  into  account 
when  computing  the  rotation.  The  calculation  of  residuals  is  made  by 
bringing  station  coordinates,  referenced  to  the  CIO  pole,  into  the 
inertial  system.  The  transformation  between  the  earth  fixed  system 
using  the  conventional  International  Origin  (CIO)  and  the  inertial  system 
is  given  by 

(35)  X £ pr  = ABCDX  inertial 

where  D = general  precession 
C = nutation 

B = rotation  from  true  inertial  equator  and  equinox  of  a 
given  time  to  Greenwich  at  that  time 
A = polar  motion  to  the  CIO  pole  using  polar  motion  values 
routinely  solved  for  in  Defense  Mapping  Agency  Navigation 
Satellite  processing. 

FORCE  EQUATION 


The  force  equation  is  represented  as 

(36)  x = G(x,*,t)  = Ae  + As  + Am  + Ad  + Ar  + A^s  + Atm  + At 

where  the  accelerations  are  due  to  earth,  sun  and  moon  gravity, 
atmospheric  drag,  radiation  pressure,  solar  and  lunar  tidal  distortion, 
and  vehicle  thrust. 

A (EARTH  GRAVITY) 


The  earth's  potential  in  an  earth  fixed  frame  is 

(37)  V = u£  2 ae  cnmPn  (rfr) cos(mx ) a§  SnmP™  (-$-)  sin(nu)! 

n=0  m=0  Vl  + )J±±L  \ 


where  r - f'  (earth  fixed  coordinates) 

w 

a„  = Semi-major  axis  of  the  earth 
= 6378.145  km 
Pm  = Legendre  polynomial 
u = Earth's  gravity  constant 
= 398601.0 

> = Longitude  with  respect  to  Greenwich 
CnrT1,  Snm  = Gravity  constants 
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Since  the  coordinate  system  has  its  center  at  the  earth's  center  of 
gravity  we  have  C1Q=C11=S11=0.  The  earth's  gravitational  acceleration 
can  now  be  given  by 


Ae  = VtV  = — 

c 1 ^\y 


where  Vi  is  the  inertial  gradient.  For  the  purpose  of  calculation  we 
rewrite  (37)  as 

(39)  V = Cm^n'CmC 

n^O  m=0 


noting  the  V°=C]i=CiQ=S]i=0. 


Define  the  transformation  E by 


(40)  E = 

3 n 


Introduce  the  longitude  by 


c(x)  = sin  (o)  cos  (a) 


/ - out  \ni  -pj-  ^_i  cli  xi 

s ( A ) = sin  (e)  sin  (x)  = ^ y E x- 

M i = ! 


where 


e = cos' 


1 /_al 


= inertial  components  of  r . 


Using  the  recurrence  relations  for  Legendre  polynomials  we  have 
recurrence  relations  for  U and  V given  by 


Cl  = 


(n-m+1)  [r 


(2  n+1 ) Up  - (n+m)  p Un_-| 


(n-m+1 ) | r 


-9i  (2  n+1  ) V™  - (n+m)  p V^_-, 


J 


F 


r 


= (2  n+1)  p[UR  C(x ) - VR  S ( X ) ] 
C]  = (2  "+1>  pfvn  C(a)  ' un  S(A)] 


Equation  (42)  is  called  Horizontal  stepping  and  (43)  is  called  Diagonal 
stepping. 


Using  the  values 

U°  = lJ/|r| 
o 

v0  = 0 


yo  = uaeg 1 
1 fr] 

V?  = 0 


we  start  at  n=l  , m=0  in  the  horizontal  stepping  equation  and  comoute  1 1 
UP,  vP  for  i=2,3...N.  We  then  utilize  diagonal  stepping  and  calculate  Ug  ,V-| 
1 1 Returning  to  norizontal  stepping  enables  the  calculation  of 

u],  vj  for  i=2,3,...N.  This  process  is  repeated  until  m = M,  where 
M N . Note  that 

N-m  _ ( - 1 )m( n-m) ! J* 

un  - un 

(n+m) ! 

v-m  = (-Dm+1(n-m)! 

n (n+m)! 


We  can  now  compute  (33)  by 
N n 

(44)  vTV  = E*  2 L [CnmaeVU™  + snmaeWR] 
n=0  m=0 


where  v is  the  earth  fixed  gradient. 

The  recurrence  relations  for  U and  V can  be  given  by 
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(45) 


m m- 1 m+1 

% An  Dp+1  Un+] 

m m-1  m+1 

-%A n Vn+1  -%  Vn+1 


-(n-m+1 ) U 


n+1  J 


aev^n 


v yin-1  _ ! vm+l 


2 n n+1 


n+1 


% A™  U^j  +%  u'nntj 


(n-m+1 )Vm 


n+1 


where  Ap  = (n-m+1)  (n-m+2) 


A , A (SUN  AND  MOON  GRAVITY) 

Coordinates  of  the  sun  and  moon  are  stored  on  the  Sun-Moon  file 
at  one  day  and  one  half  day  intervals  respectively.  A sixth  order 
Lagrangian  interpolation  procedure  is  used  to  obtain  the  values  at 
any  time.  The  sun's  acceleration  on  the  satellite  is  given  by 

(46)  As  = -y  s f r~rS  + — — ^ , ys  = .1330614  E(12) 

\l  r-rs | 3 |rs|3y 

The  expression  for  the  moon  is  similar  with  ym  = .490074  E( 04 ) . 


Ats’  Atm  (TIDAL  DISTORTION) 

The  gravitational  (tidal)  attraction  of  the  sun  and  moon  causes  the 
earth  to  become  elongated  on  an  axis  pointing  toward  the  disturbing  body. 
The  redistribution  of  mass  results  in  a perturbation  of  the  earth's  own 
gravitational  field,  which  is  represented  by  the  potentials 


(47) 

where  k|_ 


= -k|  rV 

L I rs  I 3 


aP5 


P2(r* 


-k, 


ym 


aP5 


. rm  i 3 1 rT3^ 


is  Love's  constant  and  P2  is  the  Legendre  polynomial. 


^2 ( r*Hn) 


The  associated  acceleration  on  the  satellite  is  given  by 
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(48)  A 


Ms 


ts 


,l3 


ae5 

Ms 


C(-1l(rVs)2+3/2)r+3(r*rJrs] 

2 


with  a similar  expression  for  Atm. 
Ar  (RADIA1 ION  PRESSURE) 


A shadow  test  is  performed  to  determine  if  the  satellite  is  in 
sunlight  or  shadow.  If  r-rs>0  then  the  satellite  is  in  sunlight.  If 
r-rs<0  then  we  compute  |rxrs|.  If  |rxrs  is  less  than  ae  the  satellite 
is  in  shadow  and  if  not  then  it  is  in  sunlight.  When  the  satellite 
is  in  sunlight  we  compute  the  radiation  pressure  acceleration  by 

(49)  Ar  = lv|  1014  (r-rs) 

[T-rsP 

where  s is  the  satellite  cross-sectional  area  and  m is  satellite  mass. 
Ar  is  set  to  zero  in  the  shadow. 

Ad  (ATMOSPHERIC  DRAG) 


The  relative  velocity  of  the  satellite  with  respect  to  the 
atmosphere  is 


Vr  = r-uxr 


where  [ r 
(CD)*/ 


,r)  is  the  inertial  satellite  position,  velocity  and  cr  is 
. The  acceleration  of  the  satellite  due  to  drag  is 


(50)  Ad  = -CDP  — I Vr|  Vr 
2m 

where  m is  the  satellite  mass,  s is  the  cross-sectional  area  and  p the 
atmospheric  density  function  defined  by 

P = Exp  [Ah-B- (Ch2+Dh-E ) \ 

The  height  h is  the  satellites  geocentric  altitude  above  its  subpoint 
given  by 


- - 1 -• 
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1 


h 


r|  -R 


where  R is  the  distance  from  the  earth's  center  to  the  satellites  sub- 
point. 


R 

e 


ae 


/l  + sin^  ej^ 
geocentric  latitude 


At  (VEHICLE  THRUST) 

(51)  At  = RA 
where 

R = [r,r,rxr] 

A = 

The  values  for  A are  given  by  input  as  constants  in  the  radial, 
tangential  and  out  of  plane  directions. 


VARIATIONAL  EQUATIONS 


Partial  derivatives  are  computed  by  the  same  integration  procedure 
simultaneously  with  the  orbit  integration.  The  variational  equations 
are 


(52) 

where 


h = 


3G 

3x 


+ 


_G 

9p 


0 p = osculating  orbital  element 
aG  . J 3 Model  p=CD,  A or  kr 

3p  ] 3 Model  Parameter 


V, 


24 


EARTH  GRAVITY  VARIATION 


L 


(53) 

where 


''A, 


>jX 


e = v^V  = 


N 


E*  V j [C 


rf=t)  m=0 


nm 


a£  v2UlJ)  + 


2 2 m 
ae*  Un 


% A^vu"1;] 


-%  a„vU 


m+ 1 


n+1 


1 ni.m-l 
-2  Anae/Vn+1 


-%  aeW^j 


- (n-m+1 )aevll 


n+1 


(54) 


a2v2Vm  = 
de  vn 


% A”  aeW-J  aevv'n»;l 


-+  VCl  ^ aevU»*J 


-(n-m+1 )a0vV„. 


e n+1 


The  partial s for  Cnm  and  Snm  are 


3VjV 

L °e“un 


* E*aPvll?1 


(n,m)/(0,0) 


(55) 


3V,  V 

— — = E*aQWm 

c OgHn 


3Snm 


3Vj  V 


= VtV 


■7„ 


3u 
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Snma|v^Vp]E 
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SUN  AND  MOON  GRAVITY  VARIATIONS 


RADIATION  PRESSURE  VARIATIONS 


Setting  p = | r-rs | the  derivatives  are 


krS  10 


VA  = krs  10 
VHr  ,y  z~ 

m 

krS  1014 


-3(x-xs)2+  p2  \ 
-3(x-xs)  (y-ys)  j 
-3(x-xs)  (z-zs)/ 

-3(y-ys)  (x-xs A 
-3(y-ys)2  +P l 
-3(y-ys)  (z-z s)/ 

-3(z-zs)  (x-x  )\ 
-3(Z-2SS)  (y-ys5) 
-3(z-zs)2+  p2s  / 


aAr  = S 1014  (r-r<- ) 


DRAG  VARIATIONS 


TpIVJ  + pvr  3 1 vr  I + I Vr  | Vrrr^v, 
a(rr)  2m  L 3(r*)  H7t7  {rn 


3Ad  = -_£S  |Vr|Vf 
2m 


Since  Vr  = r - u.«xr=r-  Or 


n = (CD)*  w 0 0 

0 0 0 


u 0 

0 o (CD) 


•r  = I 


d^rm  * [-0,1] 

9(rr) 


1 


d | Vr | aVr 

= v*  

r a(rr) 


3(rf) 

IP 


3(rr) 

9p 


3r 

3p 

3h 

3h 

3r 


r 3p 

3r 

3p  3h 

3h  3r 


=P[A- (ch  + D/2)  (ch^+Dh-E)”^] 
I r I - h ) 3 ^ 


= 1 


[h  + 


*§ 


]r*  + 


h)3  e2  (0,0, Z) 


| r | 2 a?  (1-e2) 


THRUST  VARIATIONS 


(60) 


3A+ 


3R 


3 ( rr ) 
3A 


3 (rr) 


Setting 


3R 

s(rf') 


3r 


3f 


3rxf 


3r 

aTrfT 


3(rf')  ’3(rf)  *9(rf) 

r 


3£ 

3r 


3r 

3r 


3r 


3 r 
3r 


= 1 


_3r 

[I  - r r*] 


, 0 


3T 


WF)  = Co’  !£] 
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A = 

1 -u)3 

u3  1 

C\J  f— 

3 3 

l 

u>i  - Aq 

<x>2 

= Ap 

-Uj2  ^1 

1 J 

0)3  = w 

(At 

+ tAt) 

t = time  in  sec. 

from  the  beginning 

of 

the  year. 

8rs  = (BCD)* 
ap 


s3 

|-sl 


-s3 

0 


u>S2 

-Us  i 
0 


wtS2 
-aits  1 
0 


=(BCD)*Q 


where  P = (Aq,  Ap,  At,  At)  is  the  polar  motion  set. 


3Dt  3D.^  9rs 

3p  3rs  9p 

3Dt 

= (BCD)*A*AQ 

3rS 

3Dt  3rs 
= AQ 


r)  D 4- 

(61)  — = — - AQ 
3P  3rSo 

Thus 

3D£  3D.  3D*  3D. 

Bpp(t)  = — = Q*(t)  A*(t)  -1  — l-  A(t)  Q(t) 

3p  3p  3rSo  3rSo 

= Cjtt)A*(t)Bss(t)  A ( t ) Q(t) 

Evaluating  AQ  at  TCA  of  a pass  gives 

(62)  Bpp  = Spass  BPP(t)  = Q*(TCA)  A*(TCA ) Bss  A (TCA ) Q(TCA) 


Equation  (62)  and  other  similar  relations  are  now  used  to  convert  the 
station  coordinate  section  of  the  pass  normal  equations  into  a polar 
motion  section.  This  technique  is  used  to  incorporate  polar  motion 
parameters  into  Celest. 
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